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Abstract 

A perturbation method for computing quick estimates of the echo decay in pulsed spin echo gradient NMR diffusion experiments 
in the short gradient pulse limit is presented. The perturbation basis involves (relatively few) dipole distributions on the boundaries 
generating a small perturbation matrix in 0(s 2 ) time, where s denotes the number of boundary elements. Several approximate 
eigenvalues and eigenfunctions to the diffusion operator are retrieved. The method is applied to 1-D and 2-D systems with Neumann 
boundary conditions. 
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perturbation method, geometry, Poisson's equation 
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1. Introduction 



NMR-methods provide an arsenal of tools to study restricted 



diffusion 
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where not only mass transportal proper- 
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ties such as flow and diffusion can be studied |4 |5|, |a |7|] 
but also characteristics of the material JH [Tol IT]]. Com- 
monly used for diffusion studies is the pulsed gradient spin- 
echo (PGSE) NMR technique where the particle positions are 
labeled by a magnetic field gradient |12]. Position labeling 
is commonly performed by finite-length magnetic field gra- 
dient pulses and the theory for this experiment is described 
by the Bloch-Torrey equations lfl3ll . In the short time gradi- 
ent limit however, the spin-echo decay simplifies to a Fourier 
transform over the propagator [3]. The short time gradient ap- 
proximation (SGP) is commonly used to describe the diffusion 
process when the geometric length scales of the material are 
longer that the effective gradient length scales as given in the 
q-vector approach lllil Il5[ [la, flTII . In heterogeneous materials 
the spin-echo decay normally results in a function that can be 
describe by a sum of exponentially decaying functions result- 
ing in a rather featureless form. However, in structurally well- 
defined materials, such as packed mono-disperse micrometer 
sized beads it can display detailed features from which mate- 
rial structure details can be obtained [0 001 . 
In addition, in the limit of short gradient pulses, the initial 
slope of the spin-echo decay always conveys information of the 
mean square displacement independent on material homogene- 
ity/heterogeneity. It is thus of interest to calculate a spin-echo 
decay from homogeneous and heterogeneous materials in order 
to learn more about the dependence between material structure 
and diffusion. The naive approach to calculate the echo decay 
in the SGP-limit is done by diagonalizing the diffusion operator. 



In this paper we develop a perturbation technique to calculate 
rough, but quick estimates of the echo decay, based on approx- 
imate eigenfunctions of the diffusion operator. These approx- 
imate eigenfunctions separates free diffusion and the influence 
of the material. Interesting features of the material can thus be 
analyzed in detail, also for large scale systems. 

2. Theory 

In the short gradient pulse approximation the echo decay is 
described by the so called master equation, a Fourier transform 
over the propagator 12, 221 



E(q,t) = -(q\P(r,r',t)\q) 



(1) 



where by the volume term i we have assumed that the initial 
positions of the particles are equally distributed among the sam- 
ple and (q\ = e' 2nqz , where q is a real and the applied gradient 
is in the z-direction. The propagator in equation Q] denotes the 
ordinary diffusion propagator 12211 . which can be expanded in 
eigenfunction/eigenvalue pairs |23j,|22[] as 



i=0 
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The eigenequation for the eigenfunction/eigenvalue pairs is 
written as 

L\i) = At\i) (2) 

where L denotes the effective diffusion operator associated with 
the boundary conditions, which will be defined in detail below. 
Note that the master equation (equation [TJ can be written as 
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i.e. the echo decay is defined by the overlap between incoming 
modes (q\ and the eigenfunctions |/) to the diffusion operator L, 
weighted by the time-dependent term e~ a> . In general, equation 
[3] must be solved using numerical methods, since the eigen- 
functions of the diffusion operator are known only for simple 
geometries. We note also that for periodic boundary conditions 
the incoming mode (q\ is an harmonic function e.g an eigen- 
function to A when q is an integer. 
We write the diffusion operator as 

L = A - S (4) 

where A denotes the ordinary Laplace operator and S denotes 
an operator defining the boundary conditions. We will refer to 
S as the surface operator. We note that the unperturbed prob- 
lem reduces to Laplace equation with solutions (q\ if q is a valid 
wave number [124] and it is evident that part of the perturbation 
basis need to consist of a set of integer valued (q\, for a cor- 
rect solution in absence of S . By the form of equation [3] we 
are motivated to find some set of functions orthogonal to (q\ 
describing the influence of the surface operator S . The eigen- 
functions of the diffusion operator would then be approximated 
by linear combinations of (q\ and these unknown vectors. The 
echo decay in equation[3]would then reduce to 

E( qi ,t)*^e- U '\a qi \ 2 , (5) 

for weights a qi induced by the surface operator S . For sim- 
plicity we will omit the index i. A' denotes the approximate 
eigenvalues of L. Such a construction will now be explained in 
detail. 

We construct S by assuming Neumann conditions at the 
boundary Q 

ft- V0(weQ,O =0, (6) 

for the (unknown) solution (p(r, t). The operator S equals ft ■ V, 
and acts as a directional derivative on Q. Each eigenfunction of 
S consist of two (5(r-<x>)-functions with sign change over Q and 
it is clear that standard perturbation techniques will not work, 
as the norm of S is large in the Laplace basis. By the form of 
the eigenfunctions to S we will refer to them as dipoles. Now 
we Fourier expand the eigenspace of S , with sign change over 
Q to preserve the dipole form and denote such surface Fourier 
modes by cr s . If the surface is smooth, a Fourier expansion on 
the boundary captures incoming waves (q\ of about the same 
wave numbers. This means that for a truncated set of Fourier 
modes {(t7|}^ =1 in the low ^-regime, a set of low wave-number 
surface modes suffices. We denote the number of such surface 
modes by M. The corresponding surface functions \s), are then 
calculated by solving Poisson's equation 

\s) = / 1 o- s {u)d(D (7) 
n 

where in two dimensions the kernel is replaced by log \r - u>\. 
The approximate solutions to the diffusion problem (equation 
|2]l can then be written as linear combinations of eigenfunctions 



to the Laplace operator A, \q) and solutions \s) to Poisson's 
equation (equation|7]i 

N M 

I0 = 2>*l?> + 5>l'>- (8) 

q s 

This linear combination does not bear sense if N -»■ oo, as of 
course {(gD^o already form complete set. If we however re- 
strict ourselves to a subset of eigenfunctions of A, N < oo, the 
complementary basis spanned by \s) is interesting and proposes 
a perturbation techniqu^H The outline of the mixed basis ap- 
proach can also be found in lE5ll . 

Although the surface distributions ct s ((l>) are chosen to be 
orthogonal, the corresponding \s) will not be, but more impor- 
tantly they nor will be orthogonal to (q\. The orthogonalization 
is however straight forward noting that 

(<?k> = -^-(q\o- s ), 

since A is self-adjoint. The scalar product between two solu- 
tions to Poisson equation is also involved in the orthogonaliza- 
tion, but can be treated as follows 

J J V-u>\ J \r-co'\ 

= J cr s (a))cr s '(cD')®(tL>, a)')da)da>' . 

The interchange of the integration variables is valid provided 
that the unit cell is neutral [26] and the resulting function 
®(ai, a)') can be pre-calculated! Importantly the function 
only depends on the unit cell size and number of dimensions 
and thus needs only to be calculated once and used as a look-up 
table. Furthermore, by noting that both 

(A-S)\s ± ) 

and 

(A-S)\q) 

only need to be calculated over the surface, the resulting per- 
turbation matrix is quickly accessible and is of size (Af + M) x 
(Af + M). The perturbation matrix has the following form 

A u =(i\(A-S)\j), 

where i and j range over all perturbation vectors i.e. {(#|}^ =1 
and The scalars a q in equation [5] as well as the ap- 

proximate eigenvalues can be retrieved by diagonalization of 
the small perturbation matrix A. We denote the eigenfunctions 
of the perturbation matrix by \k). For t -* 0, the echo decay is 
solely expressed by |(t7|A;)| 2 and can directly be read out from 
the first Af elements of the diagonal of the perturbation matrix. 
Note also that when t -> oo the echo decay reduces to the first 
column of the perturbation matrix. 



' This restriction also connects naturally with the experimental NMR setup, 
where the range of (/-vectors is not complete but restricted by the gradient 
strength available. 
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3. Results 



The perturbation basis has been validated in several trivial 
and non-trivial domains with good results. Three examples are 
presented here and in all examples the free space diffusion con- 
stant is set to unity. 

The first example consists of diffusion between two plates 
separated by a distance a, a well studied situation for which an 
analytic expression is known ESQ. A standard finite dif- 
ference approach is used with a grid spacing h = a/50. The per- 
turbation basis consist of N = 10 eigenfunctions to the Laplace 
operator A, and one dipole function representing the boundary 
(M = 1). The echo decay is shown in figure[T]for t = 100 and 
t = oo together with the real SGP-signal calculated from the 
eigenfunction expansion of L and the analytical infinite time 
solution E(q,t -*■ oo) = |sinc(7rq , /)| 2 J2l- The relative error of 
the approximate echo decay is of order ~ 10~ 4 for this perturba- 
tion basis and relative error of the apparent diffusion constant, 
estimated from the initial slope of the echo decay 112911 . is also 
of order ~ 10~ 4 . 

The following two examples consist of two-dimensional sys- 
tems using 4 • 10 4 grid points. Periodic boundary conditions 
are used on the computational cell, which has a side length 
/ = 200. Neumann boundary conditions separate the void space 
(white regions) and the structure (grey regions) and the echo 
decays are calculated in the void space. The dipole distribu- 
tions for the boundaries are calculated by diagonalizing the fi- 
nite difference approximation on the boundary yielding Fourier 
modes spanning the surface and sign change over the domain 
preserve the dipole form. The first such example consist of ran- 
domly distributed discs with equal radius (see figure |2)- Fig- 
ure [3] shows the real and approximate echo decays for times 
t = 200,900,2000 which are chosen to double the diffusive 
length in each step. The real echo decay is calculated with equa- 
tion[3]using the 280 first eigenfunctions/eigenvalues to L, which 
gives an error < 10~ 9 for t > 200. The approximate echo decay 
is calculated using the first N = 150 eigenfunctions to A and 
20 surface Fourier modes per disc are used, in total M = 280 
surface functions represent the boundaries. The relative error 
of the echo decay is of order ~ 10~ 3 . 

The last example consists of diffusion in a more interest- 
ing 2-dimensional model. The model is generated by a par- 
ent/child process 113011 . where parents are created randomly us- 
ing a uniform distribution and children are distributed around 
each parent using a Gaussian distribution. The pixel positions 
of the children then represent the material. The number of 
parents/children and the distribution parameters can be varied 
and the space of geometries is rich. Although such geometries 
work well in discrete case, they can of course not be spanned 
by Fourier modes in the continuous limit. Figure |4] shows an 
example of one such geometry and echo decays calculated for 
times t = 200, 900, 2000. The perturbation matrix is calculated 
using the N = 100 first eigenfunctions to A and M = 100 surface 
vectors. In the calculation of the real echo decay the first 230 
eigenfunctions corresponding to the void space are used, which 
gives an error of order < 10~ 9 for t > 200. The relative error of 
the approximative echo decay is of order 10~2. 
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Figure 1 : The figure shows echo decays for diffusion between 
two plates, separated by distance a. The real echo decay is 
calculated using equation [3] with the full spectrum of the dif- 
fusion operator L for time t = 100 (•). The approximate echo 
decay for the corresponding time is calculated using N = 10 
eigenfunctions to the Laplace operator A and M = 1 surface 
function (solid line). Also shown is the infinite time solution 
E(q, oo) = |sinc(7rg/)| 2 (■) and the approximate infinite time 
solution (dashed line) using the same perturbation basis as for 
the t = 100 signal. The approximate signals coincide well with 
expected results (the relative error is of order 10" 4 ). 




Figure 2: The figure shows a 2D-system consisting of randomly 
distributed discs of equal radius. The system consist of 4 • 10 4 
grid points and Neumann boundary conditions separates the 
void space (white region) from structure (grey region). Figure|3] 
shows the real and approximate echo decay for different times. 
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Figure 4: Echo decays for the model shown in top left image for times t = 200 (top right), t = 900 (bottom left) and t = 2000 
(bottom right). The times are chosen to approximately double the diffusive length for each step. The real echo decay E(q,t) 
is calculated using equation [3] (•) and the approximative echo decay is calculated using M = 100 surface functions and the first 
N = 100 eigenfunctions of A (solid line). Note that the approximative echo is calculated only at the corresponding ^-values but the 
solid line is drawn between these, for visualization. I = 200 denotes the box side length. The relative error of the approximative 
echo decay is of order ~ 10" 2 . 
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4. Conclusions 




Figure 3: The figure shows echo decays for the 2D-example of 
randomly distributed discs (see figure |2}- The real echo decay 
(calculated using equation [3} is shown for times t = 200 (▼), 
t = 900 (■) and t = 2000 (A). The approximate echo decay is 
calculated using N = 150 eigenfunctions to the Laplace opera- 
tor A and M = 280 surface functions and is shown for t = 200 
(dashed line), t = 900 (dotted line) and t = 2000 (filled line). 
The box side length is I = 200. The times are chosen to approx- 
imately double the diffusive length for each time and the rela- 
tive error of the approximate echo decays is of order ~ 10~ 3 . 
Note that the approximative echos are calculated only at the 
corresponding ^-values but lines are drawn between these, for 
visualization. 



We have shown that for diffusion problems with Neumann 
boundary conditions the echo decay in the SGP-limit can be cal- 
culated via a perturbation method with a mixed basis. Approx- 
imate echo decays are presented together with analytic and real 
echo decays (calculated from the eigenfunction expansion of 
the diffusion operator) for trivial and non-trivial geometries and 
the relative error of the echo decay is small. The mixed basis 
consist of (analytically known) eigenfunctions to the Laplace 
operator and solutions to Poisson's equation with dipole distri- 
butions on the boundary. Relatively few base vectors are needed 
for good result, resulting in a quickly accessible perturbation 
matrix. The method is formulated on the boundary, apart from 
a volume dependent function, which however is geometry inde- 
pendent and can be pre-calculated using standard Ewald sum- 
mation techniques, saved to disk, and used as a lookup table 
for arbitrary geometries. This reduces the calculations of ap- 
proximate propagators and/or echo decays in the SGP-limit to 
a computational complexity of 0(s 2 ), where s is the number of 
surface elements. 

As the approximate eigenfunctions not fully compensate for 
the Neumann conditions on the boundaries a resonance effect 
has been observed when using harmonic functions in the per- 
turbation basis with wave-lengths corresponding to the struc- 
ture domains (grey regions), this increases the error of the echo 
decay at ^-values corresponding to such wave lengths. At such 
wave lengths the approximate eigenfunctions consist of lin- 
ear combinations of eigenfunctions corresponding to the outer 
(white region) and inner (grey regions) domains. This effect 
can be minimized by increasing the number of surface modes 
M, preserving the orthogonality to the inner domains and or not 
introduce harmonic functions at the resonance points. Note that 
the error due to this resonance effect is of the same order as the 
error at other g-values in the geometries presented, but more 
pronounced. 

The method share similarities with other methods formulated 



on the boundary such as the boundary element methods Oi l 
(BEM), analytic element methods 0211 (AEM) and boundary 
approximation methods J33, 34] (BAE), also known as Trefftz 
methods, but might be an alternative due to the small size of the 
resulting perturbation matrix achieved in 0(s 2 ) time where s 
is the number of boundary elements. The approximate signals 
can also be improved by using the approximative eigenfunc- 
tions/eigenvalues as an initiator for other iterative method as for 
example [35] which relies on an initial guess of the eigenval- 
ues. The mixed basis approach can also be extended from the 
SGP-limit to cover time-dependent gradients using the matrix 
formulation developed by Callaghan 113611 based on a multiple 
propagator approach Q~7ll . As the standard methods for calcula- 
tions of the diffusion propagator are impractical for large-scale 
systems, due to the heavy computational demand, the mixed 
basis approach is suggested as a realistic tool for calculating 
approximative echo decays, also for finite gradients. 
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